clear; clc; close all;

%% 1. 设置参数
% 初始条件
X0 = [-2, -0.2, -1, -1, -1.12, 0.1, 0.3];

% HR 神经元1 参数
params.a1 = 1.5;    % 激活参数：影响膜电位与恢复变量之间的相互作用
params.b1 = 3;      % 非线性抑制系数：调节膜电位的非线性响应
params.c1 = 1;      % 恢复变量：决定恢复变量的响应特性
params.d1 = 5;      % 抑制作用：对 HR 神经元1的膜电位平方项产生抑制
params.i1 = 0.45;   % 外部输入电流/偏置

% FN 神经元 参数（FitzHugh-Nagumo 模型）
params.a2 = 0.77; 
params.b2 = 0.3;
params.c2 = 0.8;
params.i2 = 0;

% HR 神经元2 参数
params.a3 = 0.9;
params.b3 = 3;
params.c3 = 1;
params.d3 = 5;
params.i3 = 0.5;

% 耦合系数
params.m12 = 0.785;
params.m21 = 0.52;
params.m23 = 0.785;
params.m32 = 1.01;

% 固定参数：FN 神经元参数 ε
eps_fixed = 4.87;

% 定义要扫描的参数范围
eps_range = linspace(1, 8, 100);

% 积分与 LE 计算参数
T_trans = 300;
T_le = 500;
dt_save = 20;
RelTol = 1e-6;
AbsTol = 1e-8;

numLE = 2;  % 计算前两个 Lyapunov 指数

%% 2. 预分配存储 LE 的数组
LE1_array = zeros(length(eps_range), 1);
LE2_array = zeros(length(eps_range), 1);

%% 3. 逐点计算 LE1, LE2
for idx = 1:length(eps_range)
    eps_val = eps_range(idx);
    LE_vals = LEs(eps_val, params, X0, T_trans, T_le, dt_save, RelTol, AbsTol, numLE);
    LE1_array(idx) = LE_vals(1);
    LE2_array(idx) = LE_vals(2);
    fprintf('eps = %.3f, LE1 = %.4f, LE2 = %.4f\n', eps_val, LE_vals(1), LE_vals(2));
end

%% 4. 绘图
figure; hold on;
plot(eps_range, LE1_array, 'r-', 'LineWidth',1.5);
plot(eps_range, LE2_array, 'b-', 'LineWidth',1.5);
xlabel('\epsilon', 'FontSize',12);
ylabel('李指数', 'FontSize',12);
title('LE1 (红) and LE2 (蓝) vs. \epsilon', 'FontSize',14);
legend('LE1','LE2','Location','best');
grid on;
set(gca,'FontSize',12);